PRO MODELING_DRIVE_MORPH, havesave, drivemorph, suff, ALPHA=alpha
  
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star

  IF n_elements(havesave) EQ 0   THEN havesave   = 0b
  IF n_elements(drivemorph) EQ 0 THEN drivemorph = 0b
  IF n_elements(suff) EQ 0       THEN suff       = ''
  IF n_elements(rb3) EQ 0 THEN $
     restore,'irdc_dist_model/bgps_rb3'+suff+'.sav',/ver
  
  ;; If you're lazy and didn't notice that rb3 is new, but left
  ;;    havesave=1b, then this routine is smarter than you and runs the
  ;;    fill_in routine(s) for you. 
  IF total(rb3.morph_d) EQ 0 THEN BEGIN
     havesave = 0b
     drivemorph = 1b
  ENDIF
  
  IF ~ havesave THEN fill_in_rb3_morph,SUFF=suff,DRIVEMORPH=drivemorph,$
                                       ALPHA=alpha
  
  ;; Drive the morphological matching for the current list of sources...
  cind = WHERE(rb3.constrain, nc)
  print,nc
  ;; rb3 = rb3[cind]
  n = n_elements(rb3)
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF ~ exists THEN galactic_params 
  defsysv, '!IRDC', exists = exists
  IF ~ exists THEN irdc_params
  
  
  ;;====================
  ;; Escape sequence
  IF n_elements(rb3) EQ 1 THEN goto,eoc
  
  ;; Prep for the plotting...
  
  restore,'irdc_dist_model/morph_comparison'+suff+'.sav',/ver
  
  fitdist = rb3.fitdist
  morph_d = rb3.morph_d
  close   = rb3.equal_prob
  
  match = abs(morph_d-fitdist) LT 1.
  mismatch = ~ match
  
   ;;=========================================================================
  ;; FIGURE #0, Distance, Model vs. Morph
  myps,'./irdc_dist_model/analysis_plots/distance_model_morph.eps',$
       xsize=5,ysize=5
  
  bind = WHERE( abs(morph_d-fitdist) GE 1., nb)
  print,nb

  plotsym,0,0.5,/fill
  plot,fitdist,morph_d,psym=8,/isotropic,xtit='IRDC Modeling Distance [kpc]',$
       ytit='Morphological Matching Distance [kpc]',/nodata
  
  plots,0,1
  plots,19,20,/cont,linestyle=2,thick=2,color=cgColor('BLK5')
  plots,1,0
  plots,20,19,/cont,linestyle=2,thick=2,color=cgColor('BLK5')
  
  oploterror,fitdist[bind],morph_d[bind],morph_d[bind]*0.,rb3[bind].morph_sigm,$
             /lobar,psym=3,errcolor=cgColor('BLK4')
  oploterror,fitdist[bind],morph_d[bind],morph_d[bind]*0.,rb3[bind].morph_sigp,$
             /hibar,psym=3,errcolor=cgColor('BLK4')
  
  
  oplot,fitdist,morph_d,psym=8
  flags = WHERE(rb3.morph_constrain EQ 0,nclose)
  oplot,fitdist[flags],morph_d[flags],psym=8,color=cgColor('Red')
  
  al_legend,/top,/left,['N!dbad!n = '+string(nb,format="(I0)"),$
                        'N!dred!n = '+string(nclose,format="(I0)")],$
            box=0,charsize=0.8,textcolor=['Black','Red']
  
  myps,/done
  
  ;;=========================================================================
  ;; FIGURE #1, Distance, Model vs. Morph
  
  myps,'./irdc_dist_model/analysis_plots/dist_fitvmorph'+suff+'.eps',$
       xsize=15,ysize=9.5
  multiplot,[3,2],gap=0.04,/doxaxis,/doyaxis,mtitsize=1.0,$
            mtitle='[Morph Match - IRDC Model]'
  
  ;;=======================================
  ;; PANEL 1 -- DIST COMP
  
  plotsym,0,0.5,/fill
  plot,fitdist,morph_d,psym=8,/isotropic,xtit='IRDC Modeling Distance [kpc]',$
       ytit='Morphological Matching Distance [kpc]',/nodata
  
  plots,0,1
  plots,19,20,/cont,linestyle=2,thick=2,color=cgColor('TG5')
  plots,1,0
  plots,20,19,/cont,linestyle=2,thick=2,color=cgColor('TG5')
  
  oplot,fitdist,morph_d,psym=8
  flags = WHERE(close EQ 1,nclose)
  oplot,fitdist[flags],morph_d[flags],psym=8,color=cgColor('Red')
  
  bind = WHERE( abs(morph_d-fitdist) GE 1., nb)
  print,nb
  al_legend,/top,/left,['N!dbad!n = '+string(nb,format="(I0)"),$
                        'N!dred!n = '+string(nclose,format="(I0)")],$
            box=0,charsize=0.8,textcolor=['Black','Red']
  
  multiplot,/doxaxis,/doyaxis
  
  ;;=======================================
  ;; PANEL 2 -- VERSUS IRDC Contrast
  
  plot,rb3.c_meas,morph_d - fitdist,psym=8,xtit='IRDC Contrast',/nodata,$
       ytit='(MORPH - MODEL) distance [kpc]',yr=[-5,15],/yst
  vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),1.0
  vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),-1.0
  oplot,rb3.c_meas,morph_d - fitdist,psym=8
  oplot,rb3[flags].c_meas,morph_d[flags]-fitdist[flags],psym=8,$
        color=cgColor('Red')
  
  ;; Kick out 'odd' source
  ind = where(morph_d - fitdist LT -5, nind)
  IF nind GT 0 THEN print,rb3[ind].cnum
  
  
  
  multiplot,/doxaxis,/doyaxis
  
  ;;=======================================
 ;; PANEL 3 -- VERSUS BGPS Flux Density
  
  plothist,alog10(rb3.s_peak/1.d3),xarr,yarr,/noplot,bin=0.1
  plot,rb3.s_peak/1.d3,morph_d - fitdist,psym=8,/xst,$
       xtit='BGPS 40" Flux Density [Jy]',$
       /nodata,ytit='(MORPH - MODEL) distance [kpc]',yr=[-5,15],/yst,/xlog,$
       xr=10.^[min(xarr)-0.1,max(xarr)+0.1]
  vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),1.0,/xlog
  vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),-1.0,/xlog
  oplot,rb3.s_peak/1.d3,morph_d - fitdist,psym=8
  oplot,rb3[flags].s_peak/1.d3,morph_d[flags]-fitdist[flags],psym=8,$
        color=cgColor('Red')
  
  multiplot,/doxaxis,/doyaxis
  
  ;;=======================================
  ;; PANEL 4 -- VERSUS d_model / d_tan
  
  tandist = !MW.R0*cos(rb3.l * !dtor)/1.d3/cos(rb3.b * !dtor)
  
  plot,morph_d / tandist,morph_d - fitdist,psym=8,$
       xtit='d!dMorph Match!n / d!dtan!n',xr=[0,2.5],$
       /nodata,ytit='(MORPH - MODEL) distance [kpc]',yr=[-5,15],/yst
  vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),1.0
  vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),-1.0
  oplot,morph_d / tandist,morph_d - fitdist,psym=8
  oplot,morph_d[flags]/tandist[flags],morph_d[flags]-fitdist[flags],psym=8,$
        color=cgColor('Red')
  
  multiplot,/doxaxis,/doyaxis

  ;;=======================================
  ;; PANEL 5 -- VERSUS vlsr
  
  vtan = vphys2vlsr(rb3.l,rb3.b,tandist*1.d3, VS=0., /CLEM)
  
  plot,rb3.vlsr/vtan,morph_d - fitdist,psym=8,$
       xtit='V!dLSR!n / V!dtan!n',$
       /nodata,ytit='(MORPH - MODEL) distance [kpc]',yr=[-5,15],/yst
  vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),1.0
  vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),-1.0
  oplot,rb3.vlsr/vtan,morph_d - fitdist,psym=8
  oplot,rb3[flags].vlsr/vtan[flags],morph_d[flags]-fitdist[flags],psym=8,$
        color=cgColor('Red')
  
  myps,/done,/mp
  

  ;;================
  ;; More fancy footwork
  fullrb3 = rb3
  fullfit = fitdist
  fullmor = morph_d
  
  rb3     = rb3[bind]
  fitdist = fitdist[bind]
  morph_d = morph_d[bind]
  save,rb3,fitdist,morph_d,filename='./irdc_dist_model/bgps_rb3_modmorbad'+$
       suff+'.sav',/ver
  
  openw,lun,'./irdc_dist_model/morph_mismatch.txt',/GET_LUN
  printf,lun,'     CNUM   FITDIST   MORPHDIST   COMMENTS'
  FOR ij=0L,n_elements(rb3)-1 DO BEGIN
     printf,lun,rb3[ij].cnum,fitdist[ij],morph_d[ij]
  ENDFOR
  close,lun
  free_lun,lun
  
  rb3     = fullrb3
  fitdist = fullfit
  morph_d = fullmor
  
  
  ;;=========================================================================
  ;; FIGURE #2, Distance, Model vs. Morph SUCCESS RATES  (C & S)
  
  myps,'./irdc_dist_model/analysis_plots/dist_morph_success'+suff+'.eps',$
       xsize=10,ysize=4.5
  multiplot,[2,1],xgap=0.06,/doxaxis,/doyaxis,$
            mtitle='[Morph Match - IRDC Model] Success Rates',mtitsize=0.8
  
   
  ;; Make bar chart -- IRDC Contrast
  ;; SR_BARS, array, minval, maxval, step, xp, yp, yp2
  minval = 0.d
  maxval = 0.5d
  step   = 0.05d
  sr_bars, rb3.c_meas, match, minval, maxval, step, $
           xp, yp, yp2
  
  plot,xp,yp,psym=10,xtit='IRDC Contrast',$
       ytit='Success Rate',yst=8,xmargin=[10,10],$
       thick=4,xr=[0,maxval],/xst
  axis,yaxis=1,yr=[0,250],ytit='N per '+string(step,format="(F0.2)")+' bin',$
       /save,/yst
  oplot,xp,yp2,psym=10,color=cgColor('BLK5'),linestyle=3
  
  al_legend,/top,/left,['Success Rate','N per bin'],linestyle=[0,3],$
            thick=[4,0],color=['Opposite','BLK5'],box=0,charsize=0.8,$
            linsize=0.6
  
  multiplot,/doxaxis,/doyaxis
  
  ;; SR_BARS, array, minval, maxval, step, xp, yp, yp2
  minval = -1.4d
  maxval = 0.5d
  step   = 0.1d
  sr_bars, alog10(rb3.s_peak/1.d3), match, minval, maxval, step, $
           xp, yp, yp2
  
  ;; PLOT
  plot,10.^xp,yp,psym=10,xtit='BGPS 40" Flux Density [Jy]',$
       ytit='Success Rate',yst=8,xmargin=[10,10],$
       thick=4,xr=10.^[minval,maxval],/xst,/xlog
  axis,yaxis=1,yr=[0,250],ytit='N per log(S)='+string(step,format="(F0.1)")+$
       ' bin',/save,/yst
  oplot,10.^xp,yp2,psym=10,color=cgColor('BLK5'),linestyle=3
  
  al_legend,/top,/left,['Success Rate','N per bin'],linestyle=[0,3],$
            thick=[4,0],color=['Opposite','BLK5'],box=0,charsize=0.8,$
            linsize=0.6
  
  myps,/done,/mp
  
  
  ;;=========================================================================
  ;; FIGURE #3, Distance, Model vs. Morph SUCCESS RATES  (I_MIR & GLON)
  
  myps,'./irdc_dist_model/analysis_plots/dist_morph_success2'+suff+'.eps',$
       xsize=10,ysize=4.5
  multiplot,[2,1],xgap=0.06,/doxaxis,/doyaxis,$
            mtitle='[Morph Match - IRDC Model] Success Rates',mtitsize=0.8
  
  ;; Make bar chart -- I_MIR
  ;; SR_BARS, array, minval, maxval, step, xp, yp, yp2
  minval = 0.d
  maxval = 200.d
  step   = 10.d
  sr_bars, rb3.i_mir, match, minval, maxval, step, $
           xp, yp, yp2
  
  plot,xp,yp,psym=10,xtit='I!dMIR!n [MJy/sr]',$
       ytit='Success Rate',yst=8,xmargin=[10,10],$
       thick=4,xr=[0,maxval],/xst
  axis,yaxis=1,yr=[0,250],ytit='N per '+string(step,format="(I0)")+$
       ' MJy/sr bin',/save,/yst
  oplot,xp,yp2,psym=10,color=cgColor('BLK5'),linestyle=3
  
  al_legend,/top,/left,['Success Rate','N per bin'],linestyle=[0,3],$
            thick=[4,0],color=['Opposite','BLK5'],box=0,charsize=0.8,$
            linsize=0.6
  
  multiplot,/doxaxis,/doyaxis
  
  ;; Make bar chart -- GLON
  ;; SR_BARS, array, minval, maxval, step, xp, yp, yp2
  minval = 20.d
  maxval = 70.d
  step   = 5.d
  sr_bars, rb3.l, match, minval, maxval, step, $
           xp, yp, yp2
  
  ;; PLOT
  plot,xp,yp,psym=10,xtit='GLON [deg]',$
       ytit='Success Rate',yst=8,xmargin=[10,10],$
       thick=4,xr=[maxval,minval],/xst,yr=[0,1]
  axis,yaxis=1,yr=[0,250],ytit='N per '+string(step,format="(F0.2)")+$
       ' deg bin',/save,/yst
  oplot,xp,yp2,psym=10,color=cgColor('BLK5'),linestyle=3
  
  al_legend,/top,/left,['Success Rate','N per bin'],linestyle=[0,3],$
            thick=[4,0],color=['Opposite','BLK5'],box=0,charsize=0.8,$
            linsize=0.6
  
  myps,/done,/mp
  
  
  ;; Send junk to Screen:
  print,''
  print,'Suffix: ',suff
  print,'ALL: '
  str = string(n_elements(rb3),total(mismatch),$
               total(match)/n_elements(rb3)*100.,$
               format="('N = ',I4,'   N_BAD = ',I4,"+$
               "'   Success Rate = ',F0.1,'%')")
  print,str
  
  print,''
  print,'C GEQ 0.1'
  junk = WHERE(rb3.c_meas GE 0.1, n1)
  str = string(n1,total(mismatch[junk]),total(match[junk])/float(n1)*100.,$
               format="('N = ',I4,'   N_BAD = ',I4,'   "+$
               "Success Rate = ',F0.1,'%')")
  print,str
  
  print,''
  print,'C GEQ 0.15'
  junk = WHERE(rb3.c_meas GE 0.15, n2)
  str = string(n2,total(mismatch[junk]),total(match[junk])/float(n2)*100.,$
               format="('N = ',I4,'   N_BAD = ',I4,'   "+$
               "Success Rate = ',F0.1,'%')")
  print,str
  
  
  
  GOTO,eoc                      ; Just bypass all the below stuff
  
;; ;;=========================================================================
  
;; myps,'./irdc_dist_model/analysis_plots/morph_vs_grs'+suff+'.eps',$
;;      xsize=10,ysize=15

;; delta = cgSymbol('Delta')
;; sun   = cgSymbol('Sun')

;; pind = WHERE(dsep GE 0., np)
;; print,np
;; dsep    = dsep[pind]
;; vsep    = vsep[pind]
;; rb3     = rb3[pind]
;; grsi    = grsi[pind]
;; fitdist = fitdist[pind]
;; morph_d = morph_d[pind]
;; grs     = grs[grsi]

;; multiplot,[2,3],gap = 0.04,/doxaxis,/doyaxis

;; ;; ;; --FIRST PANEL -- Source offset from GRS (l,b,v)

;; ;; plotsym,0,0.5,/fill
;; ;; plot,vsep,dsep,psym=8,xtit='(GRS-BGPS) '+delta+'V [km/s]',$
;; ;;      ytit='(GRS-BGPS) '+delta+'D [deg]',tit='Selection Criteria'

;; ;; multiplot,/doxaxis,/doyaxis

;; ;; FIRST PANEL -- MODEL FIT vs. GRS

;; plot,grs.dis,fitdist,psym=8,title='GRS vs. Model Fit Dist',xtit='GRS [kpc]',$
;;      ytit='IRDC model fit [kpc]',/isotropic,/nodata

;; top = !x.crange[1]

;; plots,0,1
;; plots,top-1,top,/cont,linestyle=2,thick=2,color=cgColor('TG5')
;; plots,1,0
;; plots,top,top-1,/cont,linestyle=2,thick=2,color=cgColor('TG5')

;; oplot,grs.dis,fitdist,psym=8

;; jind = WHERE( abs(grs.dis-fitdist) GT 1., nbd)
;; cgtext,0.5,top-1,'N!dbad!n = '+string(nbd,format="(I0)"),charsize=0.8,$
;;        color=cgColor('RED6')

;; multiplot,/doxaxis,/doyaxis

;; ;; SECOND PANEL -- MORPH MATCH vs. GRS

;; plot,grs.dis,morph_d,psym=8,title='GRS vs. Morph Match Dist',xtit='GRS [kpc]',$
;;      ytit='Morph Match [kpc]',/isotropic,/nodata

;; top = !x.crange[1]

;; plots,0,1
;; plots,top-1,top,/cont,linestyle=2,thick=2,color=cgColor('TG5')
;; plots,1,0
;; plots,top,top-1,/cont,linestyle=2,thick=2,color=cgColor('TG5')

;; oplot,grs.dis,morph_d,psym=8

;; jind = WHERE( abs(grs.dis-morph_d) GT 1., nbd)
;; cgtext,0.5,top-1,'N!dbad!n ='+string(nbd,format="(I0)"),charsize=0.8,$
;;        color=cgColor('RED6')

;; multiplot,/doxaxis,/doyaxis

;; ;; ;; --FOURTH PANEL -- MODEL-GRS vs. DELTA D


;; ;; plot,dsep,model_del,psym=8,xtit='(GRS-BGPS) '+delta+'D [deg]',$
;; ;;      ytit='d(Model - GRS) [kpc]',/nodata
;; ;; vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),1.0
;; ;; vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),-1.0
;; ;; oplot,dsep,model_del,psym=8


;; ;; multiplot,/doxaxis,/doyaxis

;; ;; ;; --FIFTH PANEL -- MODEL-GRS vs. DELTA V

;; ;; plot,vsep,model_del,psym=8,xtit='(GRS-BGPS) '+delta+'V [km/s]',$
;; ;;      ytit='d(Model - GRS) [kpc]',/nodata
;; ;; vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),1.0
;; ;; vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),-1.0
;; ;; oplot,vsep,model_del,psym=8

;; ;; multiplot,/doxaxis,/doyaxis

;; ;; THIRD PANEL -- MODEL-GRS vs. GLON

;; model_del = fitdist - grs.dis
;; morph_del = morph_d - grs.dis

;; plot,rb3.l,model_del,psym=8,xtit='GLON [deg]',$
;;      ytit='d!d'+sun+'!n (Model - GRS) [kpc]',$
;;      xr=[70,10],/nodata,/xst,yr=[-15,15],/yst,title='GRS vs. Model Fit Dist'
;; vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),1.0
;; vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),-1.0
;; oplot,rb3.l,model_del,psym=8

;; multiplot,/doxaxis,/doyaxis

;; ;; FOURTH PANEL -- MORPH-GRS vs. GLON

;; plot,rb3.l,morph_del,psym=8,xtit='GLON [deg]',$
;;      ytit='d!d'+sun+'!n (Morph - GRS) [kpc]',$
;;      xr=[70,10],/nodata,/xst,yr=[-15,15],/yst,title='GRS vs. Morph Match Dist'
;; vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),1.0
;; vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),-1.0
;; oplot,rb3.l,morph_del,psym=8

;; multiplot,/doxaxis,/doyaxis

;; ;; FIFTH PANEL -- MODEL-GRS vs. GLAT

;; plot,rb3.b,model_del,psym=8,xtit='GLAT [deg]',$
;;      ytit='d!d'+sun+'!n (Model - GRS) [kpc]',$
;;      /nodata,yr=[-15,15],/yst,title='GRS vs. Model Fit Dist'
;; vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),1.0
;; vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),-1.0
;; oplot,rb3.b,model_del,psym=8

;; multiplot,/doxaxis,/doyaxis

;; ;; SIXTH PANEL -- MORPH-GRS vs. GLAT

;; plot,rb3.b,morph_del,psym=8,xtit='GLAT [deg]',$
;;      ytit='d!d'+sun+'!n (Morph - GRS) [kpc]',$
;;      /nodata,yr=[-15,15],/yst,title='GRS vs. Morph Match Dist'
;; vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),1.0
;; vline,/horiz,linestyle=2,thick=2,color=cgColor('TG5'),-1.0
;; oplot,rb3.b,morph_del,psym=8

;; multiplot,/reset
;; multiplot,/default
;; multiplot,/reset


;; myps,/done

;; fullrb3  = rb3
;; fullmord = morph_d
;; fullfitd = fitdist
;; fullgrs  = grs

;; print,n_elements(rb3)
;; jind = WHERE( abs(grs.dis-fitdist) GT 1., nbd) & print,nbd


;; rb3mismatch = rb3[jind] 

;; gind = MISSING(jind,lindgen(n_elements(rb3)),ng) & print,ng

;; rb3yesmatch = rb3[gind]

;; ;; Write new RB3 files

;; save,rb3,filename='./irdc_dist_model/bgps_rb3_grs.sav',/ver

;; rb3 = rb3mismatch
;; save,rb3,filename='./irdc_dist_model/bgps_rb3_grsno.sav',/ver

;; rb3 = rb3yesmatch
;; save,rb3,filename='./irdc_dist_model/bgps_rb3_grsyes.sav',/ver

;; spawn,'cp ./irdc_dist_model/bgps_rb3_grs*.sav /d4/idl_lib/'+$
;;       'distance-omnibus/irdc_dist_model/.'


;; morph_d = morph_d[jind]
;; fitdist = fitdist[jind]
;; rb3 = rb3mismatch
;; grs = grs[jind]

;; openw,lun,'./irdc_dist_model/grs_mismatch.txt',/GET_LUN

;; printf,lun,'       CNUM   FITDIST    MORPH_D    GRS_DIST'
;; FOR ii=0L, nbd-1 DO BEGIN
;;    printf,lun,rb3[ii].cnum,fitdist[ii],morph_d[ii],grs[ii].dis
;; ENDFOR
;; close,lun
;; free_lun,lun

;; rb3     = fullrb3[gind]
;; fitdist = fullfitd[gind]
;; morph_d = fullmord[gind]
;; grs     = fullgrs[gind]


;; openw,lun,'./irdc_dist_model/grs_yesmatch.txt',/GET_LUN
  
;; printf,lun,'       CNUM   FITDIST    MORPH_D    GRS_DIST'
;; FOR ii=0L, n_elements(rb3)-1 DO BEGIN
;;    printf,lun,rb3[ii].cnum,fitdist[ii],morph_d[ii],grs[ii].dis
;; ENDFOR
;; close,lun
;; free_lun,lun
  
;; driver_morph,cnum_list=rb3.cnum,/save,/clem
  
  
;;=========================================================================
  eoc:
  
END
